docs/rxd-tutorials/calcium waves.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: Calcium waves"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "<p>In every cell type, intracellular Ca<sup>2+</sup> signaling plays a major role in a diverse set of activities, including process homeostasis and genetic expression. In neurons, Ca<sup>2+</sup> wave propagation has been implicated as a major component of intracellular signaling and is complementary to membrane electrical signaling via action potentials. It follows that the cell must contain tightly-controlled intrinsic mechanisms for regulating intracellular Ca<sup>2+</sup> content and waves; these include cytosolic buffers, storage in the ER, and sequestration into mitochondria. In this tutorial, we examine how some of these basic mechanisms interact to produce a Ca<sup>2+</sup> wave.</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Concept\n",
    "<p>Ca<sup>2+</sup> is well-buffered by several cytosolic proteins (e.g. calmodulin, calbindin, etc...), but must still traverse a significant distance within the neuron to reach its targets. This distance is long enough that pure diffusion would pose a temporal problem for cell communication. Instead, the cell uses oscillating Ca<sup>2+</sup> waves to harness the power of chemical signaling. These waves can travel up and down the Endoplasmic Reticulum (ER) to effect homeostatic changes. Understanding how these waves are modulated in real time by the cell poses a great difficulty because of complex spatiotemporal patterns and nonuniform distribution of the channels that control Ca<sup>2+</sup> flux across the ER. Inositol triphosphate (IP<sub>3</sub>) receptors IP<sub>3</sub>Rs), Sarco/endoplasmic reticulum Ca<sup>2+</sup>-ATPase (SERCA pump), and Ryanoine receptors (RyRs). For simplicity’s sake, we will not deal with RyRs in this tutorial.</p><p>\n",
    "Diffusible Ca<sup>2+</sup> is released directly from the ER via IP<sub>3</sub>Rs upon cooperative binding of free Ca<sup>2+</sup> and IP<sub>3</sub> to the receptor. When IP<sub>3</sub> levels are ideal, this results in Ca<sup>2+</sup>-induced Ca<sup>2+</sup> release (CICR), yielding a wave that can travel along the ER. As mentioned above, these waves are implicated in the gene expression and protein upregulation associated with neuronal plasticity by way of their interaction with the neuronal soma. The SERCA pumps within ER are normally responsible for maintaining appropriate intracellular Ca<sup>2+</sup> levels. Two models exist for wave propagation: continuous and saltatory. To understand the physical and functional differences between these models, it is ideal to compare them using Reaction-Diffusion (RxD) wtihin the NEURON environment.\n",
    "</p><p>\n",
    "The model presented in this tutorial generates Ca<sup>2+</sup> waves and is based loosely off of the <a href=\"//modeldb.yale.edu/168874\">model</a> we used in <a href=\"https://doi.org/10.1162/NECO_a_00712\">Neymotin et al., 2015</a>.</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>Python is a rich and user-friendly programming language that offers significant options for systems integration. For this reason, we will run the NEURON and RxD modules within the Python interpreter.\n",
    "</p><p>\n",
    "Like other RxD problems within NEURON, the first step is to import and load the appropriate libraries into the Python interpreter.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuron import h, rxd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is often convenient to load NEURON's <code>stdrun</code> library, which provides high level functions for controlling integration; we do so here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.load_file(\"stdrun.hoc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we define parameters for the section to be used. Like any imported module within Python, the imported libraries can be accessed by the <code>h</code> and <code>rxd</code> prefixes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sec = h.Section(name=\"sec\")\n",
    "sec.L = 101\n",
    "sec.diam = 1\n",
    "sec.nseg = 101"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above code defines a section of a neuron with Length = 101 \\({\\mu}m\\), diameter = 1 \\({\\mu}m\\), and 101 discrete segments. By chunking the neuron into a large number of segments, we increase the resolution of the output at the expense of processing power. For a smaller number of segments, the opposite is true. In larger models, it is useful to represent <code>nseg</code> as a variable – this way, if changes need to be made during any number of simulations, the single variable can be altered and will affect all subsequent iterations of <code>nseg</code>. In this example, we only reference a single section, so there is no need to elaborate.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Enable variable time step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we activate the variable time-step integration algorithm. This is a class of functions which are part of the <a href=\"http://computation.llnl.gov/casc/sundials/description/description.html\">SUNDIALS</a> package and allow variable-order, variable time-step, and implicit solving of ordinary differential equations (ODE). This is a particularly useful tool for solving a system with multiple species, whose concentrations may change suddenly from stable conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.CVode().active(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The argument of <code>active</code> specifies whether or not variable step integration will be used. To use fix step integration, pass the argument <code>False</code> instead."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters\n",
    "<p>The diffusion coefficients <code>caDiff</code> and <code>ip3Diff</code> here denote the effective diffusion rate after buffering by native cytosol proteins is taken into account.</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "caDiff = 0.016\n",
    "ip3Diff = 0.283\n",
    "cac_init = 1.0e-4\n",
    "ip3_init = 0.1\n",
    "gip3r = 12040\n",
    "gserca = 0.3913\n",
    "gleak = 6.020\n",
    "kserca = 0.1\n",
    "kip3 = 0.15\n",
    "kact = 0.4\n",
    "ip3rtau = 2000\n",
    "fc = 0.83\n",
    "fe = 0.17"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An <code>rxd.Region</code> is the volume in which species exist and reactions can occur."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cyt = rxd.Region(\n",
    "    h.allsec(),\n",
    "    nrn_region=\"i\",\n",
    "    geometry=rxd.FractionalVolume(fc, surface_fraction=1),\n",
    "    name=\"cyt\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p>This uses the rxd.Region class to define the properties for the region on the immediate interior of the membrane (<code>nrn_region = \"i\"</code>). Electrophysiologically, the only regions that matter for this simulation are the spaces directly inside the membranes. The “i” marks that we are concerned with the inner region only. For any region, a geometry may be specified. For this model, we make the assumption that ER is roughly evenly-mixed throughout the neuron segment. That is, in any given chunk of neuron, we assume there is a constant fraction of that chunk that is ER and a constant fraction that is cytosol; this is indicated using an <code>rxd.FractionalVolume</code>. The geometry argument tells NEURON what to expect the region to actually look like (i.e. a cylinder has a specific 3D spatial representation and therefore requires specific diffusion modeling for ions within that particular geometry).\n",
    "</p><p>\n",
    "We set <code>surface_fraction=1</code> because for the purposes of this model we suppose that at all points the plasma membrane is in contact with the cytosol (i.e. we assume the ER does not contact the plasma membrane and that all calcium influx from ion channels must first enter the cytosol before it can be sequestered into the ER).</p>\n",
    "<p><code>name</code> is an optional argument, but it can be queried later and some tools support using it to identify the region.</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define the ER similarly, except we omit the <code>surface_fraction</code> argument as it defaults to 0; i.e. that the ER does not contact the plasma membrane:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "er = rxd.Region(h.allsec(), geometry=rxd.FractionalVolume(fe), name=\"er\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we define the membrane separating the cytosol from the ER. <code>rxd.MultiCompartmentReaction</code>s crossing this membrane will be scaled by its surface area. Given that we are assuming that all chunks of neuron have ER and cytosol mixed together in the same way, it follows that the area of the membrane within a volume should be proportional to the size of that volume. We indicate this using the <code>rxd.DistributedBoundary</code> geometry. This is typically the right rule for separating two volumes defined as <code>rxd.FractionalVolume</code>s:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cyt_er_membrane = rxd.Region(\n",
    "    h.allsec(), geometry=rxd.DistributedBoundary(2), name=\"cyt_er_membrane\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the argument <code>2</code> to <code>rxd.DistributedBoundary</code> indicates that there is <code>2</code> μm<sup>2</sup> of ER boundary per cubic micron of volume. Ideally, this number would be based on empirical data, but for now we choose it arbitrarily."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Species"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that the definition of an <code>rxd.Species</code> has the general form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = rxd.Species(regions=None, d=0, name=None, charge=0, initial=None, atolscale=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where:\n",
    "<ul>\n",
    "    <li><code>regions</code> can either be a single region or a list of regions.</li>\n",
    "    <li><code>d</code> is the diffusion coefficient for the species; if none is specified, the species does not diffuse.</li>\n",
    "    <li><code>name</code> assigns a name to the Species that syncs with the rest of NEURON (NMODL, GUI tools, variables of the form e.g. <code>cai</code>, etc)</li>\n",
    "    <li><code>charge</code> refers to the charge of the Species; e.g. a calcium ion would have charge 2</li>\n",
    "    <li><code>initial</code> refers to the initial concentration of the Species. If <code>initial</code> is not specified and the Species has a name and exists on an rxd.Region with <code>nrn_region</code> set, then NEURON's default initialization rules apply; otherwise the default is 0.</li>\n",
    "    <li><code>atolscale</code> indicates a tolerance scale factor for variable step integration and is especially important for calcium and other substances with a low concentration</li>\n",
    "</ul>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because these are keyword arguments, they can be left out or even entered in any order."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We suppose calcium exists in both the ER and cytosol, and IP3 is present in the cytosol alone:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "caCYT_init = 0.1\n",
    "ca = rxd.Species([cyt, er], d=caDiff, name=\"ca\", charge=2, initial=caCYT_init)\n",
    "ip3 = rxd.Species(cyt, d=ip3Diff, initial=ip3_init)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>We apologize, but this tutorial is incomplete and under revision. Please continue to explore the other examples.</b> If you want to get a sense of where this is going, this is based on an <a href=\"/neuron/static/docs/rxd/tut_cawave.html\">older example</a> but that file has a number of issues."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}